{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    key                                : min               max\n",
      "=========================================\n",
      "('tracl', '<i4')                    1.000            71284.000\n",
      "('fldr', '<i4')                     231.000          481.000\n",
      "('tracf', '<i4')                    -1.000           282.000\n",
      "('ep', '<i4')                       32.000           282.000\n",
      "('cdpt', '<i4')                     1.000            284.000\n",
      "('nhs', '<i2')                      1.000            1.000\n",
      "('scalel', '<i2')                   -10000.000       -10000.000\n",
      "('scalco', '<i2')                   -10000.000       -10000.000\n",
      "('counit', '<i2')                   3.000            3.000\n",
      "('ns', '<u2')                       1501.000         1501.000\n",
      "('dt', '<u2')                       2000.000         2000.000\n",
      "('gain', '<i2')                     3.000            3.000\n",
      "('igc', '<i2')                      1.000            1.000\n",
      "('afilf', '<i2')                    207.000          207.000\n",
      "('afils', '<i2')                    298.000          298.000\n",
      "('hcf', '<i2')                      207.000          207.000\n",
      "('hcs', '<i2')                      298.000          298.000\n",
      "('year', '<i2')                     1998.000         1998.000\n",
      "('trace', '<f4', (1501,))           -0.303           1.000\n",
      "=========================================\n"
     ]
    }
   ],
   "source": [
    "from toolbox.processing import *\n",
    "from timeit import Timer\n",
    "#%ls /home/stewart/su/2d_land_data/2D_Land_data_2ms/\n",
    "file = \"/home/stewart/su/2d_land_data/2D_Land_data_2ms/su/Line_001.su\"\n",
    "#file = \"/home/sfletcher/Downloads/2d_land_data/2D_Land_data_2ms/Line_001.su\"\n",
    "#initialise file\n",
    "data, params = initialise(file)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#no coordinates in the headers, but we know energy point number and channel.\n",
    "#%cat /home/stewart/su/2d_land_data/2D_Land_data_2ms/Line_001.TXT\n",
    "#%cat /home/stewart/su/2d_land_data/2D_Land_data_2ms/Line_001.SPS\n",
    "#%cat /home/stewart/su/2d_land_data/2D_Land_data_2ms/Line_001.RPS\n",
    "#%cat /home/stewart/su/2d_land_data/2D_Land_data_2ms/header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import collections\n",
    "dmap = np.memmap(file, dtype=toolbox.typeSU(1501), mode='r')\n",
    "eps = np.unique(dmap['ep'])\n",
    "for ep in eps[:1]:\n",
    "    panel = dmap[dmap['ep'] == ep].copy()\n",
    "    panel = toolbox.agc(panel, None, **params)\n",
    "\n",
    "    trace_centers = np.linspace(1,284, panel.size).reshape(-1,1)\n",
    "    trace_width = 284/(panel.size*0.5)\n",
    "    x = panel['trace'].copy()\n",
    "    x += trace_centers\n",
    "    y = np.meshgrid(np.arange(1501), np.arange(284))[0]\n",
    "    \n",
    "    x = np.split(x.ravel(), 284)\n",
    "    y = np.split(y.ravel(), 284)\n",
    "    \n",
    "    bits = [zip(x[a],y[a]) for a in range(len(x))]\n",
    "    fig = pylab.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    \n",
    "    col1 = collections.LineCollection(bits)\n",
    "    col1.set_color('k')\n",
    "    ax.add_collection(col1, autolim=True)\n",
    "    ax.autoscale_view()\n",
    "    pylab.xlim([0,284])\n",
    "    pylab.ylim([0,1500])\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    pylab.tight_layout()\n",
    "    pylab.show()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import collections\n",
    "dmap = np.memmap(file, dtype=toolbox.typeSU(1501), mode='r')\n",
    "eps = np.unique(dmap['ep'])\n",
    "for ep in eps[:1]:\n",
    "    panel = dmap[dmap['ep'] == ep].copy()\n",
    "    panel = toolbox.agc(panel, None, **params)\n",
    "\n",
    "    trace_centers = np.linspace(1,284, panel.size).reshape(-1,1)\n",
    "    scalar = 284/(panel.size*0.5)\n",
    "    panel['trace'][:,-1] = np.nan\n",
    "    x = panel['trace'].ravel()\n",
    "    x[x < 0] = 0\n",
    "    y = np.meshgrid(np.arange(1501), np.arange(284))[0].ravel() \n",
    "    \n",
    "    zero_crossings = np.where(x == 0)[0]+1\n",
    "    zero_crossings = zero_crossings[np.diff(zero_crossings) == 1]\n",
    "    #zero_crossings = np.where(np.diff(np.signbit(x)))[0]+1\n",
    "    \n",
    "    x = ((panel['trace']*scalar)+trace_centers).ravel()\n",
    "\n",
    "    xverts = np.split(x, zero_crossings)\n",
    "    yverts = np.split(y, zero_crossings)\n",
    "    \n",
    "    \n",
    "    polygons = [zip(xverts[i], yverts[i]) for i in range(0, len(xverts)) if len(xverts[i]) > 2]\n",
    "    \n",
    "    xlines = np.split(x, 284)\n",
    "    ylines = np.split(y, 284)\n",
    "    lines = [zip(xlines[a],ylines[a]) for a in range(len(xlines))]  \n",
    "\n",
    "\n",
    "    fig = pylab.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    col = collections.PolyCollection(polygons)\n",
    "    col.set_color('k')\n",
    "    ax.add_collection(col, autolim=True)\n",
    "    col1 = collections.LineCollection(lines)\n",
    "    col1.set_color('k')\n",
    "    ax.add_collection(col1, autolim=True)\n",
    "    ax.autoscale_view()\n",
    "    pylab.xlim([0,284])\n",
    "    pylab.ylim([0,1500])\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    pylab.tight_layout()\n",
    "    pylab.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import collections\n",
    "%pylab tk\n",
    "def polytrace(data, **kwargs):\n",
    "    segs = []\n",
    "    segl = []\n",
    "    nt =  data.shape[-2]\n",
    "    for i in range(nt):\n",
    "        trace = data[i]\n",
    "        \n",
    "        line = list(zip(trace, np.arange(1501)))\n",
    "        segl.append(line)           \n",
    "\n",
    "        xx = trace  #np.ma.array(trace, mask=(trace <= 0))\n",
    "        yy = np.arange(1501) #np.ma.array(np.arange(0,1501), mask=(trace <= 0))\n",
    "         \n",
    "        curve = [(0, 0)]\n",
    "        curve.extend(list(zip(xx, yy)))\n",
    "        curve.extend([(0, 1501)])\n",
    "        \n",
    "        segs.append(curve)\n",
    "    #print segs[0] \n",
    "    print ''\n",
    "    return segs, segl\n",
    "\n",
    "\n",
    "#first lets do some checks.  does of energy points should equal number of records?\n",
    "print np.unique(data['ep']).size, np.unique(data['fldr']).size\n",
    "#no duplicates - that makes it easier.\n",
    "print 251*284\n",
    "#284 traces per shot, 2 aux traces . lets have a look\n",
    "dmap = np.memmap(file, dtype=toolbox.typeSU(1501), mode='r')\n",
    "eps = np.unique(dmap['ep'])\n",
    "for ep in eps[:1]:\n",
    "    panel = dmap[dmap['ep'] == ep].copy()\n",
    "    panel = toolbox.agc(panel, None, **params)\n",
    "\n",
    "    trace_centers = np.linspace(1,284, panel.size).reshape(-1,1)\n",
    "    trace_width = 284/(panel.size*0.5)\n",
    "    buf = panel['trace'].copy()\n",
    "    buf *= trace_width\n",
    "\n",
    "\n",
    "    segs, segl = polytrace(buf)\n",
    "    fig = pylab.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    offs = (10.0, 0.0)\n",
    "    offs = list(zip(np.arange(238), np.zeros(238)))\n",
    "    col = collections.PolyCollection(segs, offsets=offs)\n",
    "    col.set_color('k')\n",
    "    ax.add_collection(col, autolim=True)\n",
    "    \n",
    "    #col1 = collections.LineCollection(segl, offsets=offs)\n",
    "    #col1.set_color('k')\n",
    "    #ax.add_collection(col1, autolim=True)\n",
    "    ax.autoscale_view()\n",
    "    pylab.xlim([0,284])\n",
    "    pylab.ylim([0,1500])\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    pylab.tight_layout()\n",
    "    pylab.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import collections\n",
    "import matplotlib.pyplot as pylab\n",
    "\n",
    "#make some oscillating data\n",
    "panel = np.meshgrid(np.arange(1501), np.arange(284))[0]\n",
    "panel = np.sin(panel)\n",
    "\n",
    "#generate coordinate vectors.\n",
    "panel[:,-1] = np.nan #prevent wrapping when flatten 2d array\n",
    "x = panel.flatten()\n",
    "y = np.meshgrid(np.arange(1501), np.arange(284))[0].ravel() \n",
    "\n",
    "#find indexes of each zero crossing\n",
    "zero_crossings = np.where(np.diff(np.signbit(x)))[0]+1 \n",
    "\n",
    "#calculate scalar used to shift \"traces\" to plot corrdinates\n",
    "trace_centers = np.linspace(1,284, panel.shape[-2]).reshape(-1,1) \n",
    "gain = 0.5 #scale traces\n",
    "\n",
    "#shift traces to plotting coordinate\n",
    "x = ((panel*gain)+trace_centers).ravel()\n",
    "\n",
    "#split each vector at each zero crossing\n",
    "xverts = np.split(x, zero_crossings)\n",
    "yverts = np.split(y, zero_crossings)\n",
    "\n",
    "#we only want the vertices which outline positive values\n",
    "if x[0] > 0:\n",
    "    steps = range(0, len(xverts),2)\n",
    "else:\n",
    "    steps = range(1, len(xverts),2)\n",
    "\n",
    "#turn vectors of coordinates into lists of coordinate pairs\n",
    "polygons = [zip(xverts[i], yverts[i]) for i in steps if len(xverts[i]) > 2]\n",
    "\n",
    "#this is so we can plot the lines as well\n",
    "xlines = np.split(x, 284)\n",
    "ylines = np.split(y, 284)\n",
    "lines = [zip(xlines[a],ylines[a]) for a in range(len(xlines))]  \n",
    "\n",
    "#and plot\n",
    "fig = pylab.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "col = collections.PolyCollection(polygons)\n",
    "col.set_color('k')\n",
    "ax.add_collection(col, autolim=True)\n",
    "col1 = collections.LineCollection(lines)\n",
    "col1.set_color('k')\n",
    "ax.add_collection(col1, autolim=True)\n",
    "ax.autoscale_view()\n",
    "pylab.xlim([0,284])\n",
    "pylab.ylim([0,1500])\n",
    "ax.set_ylim(ax.get_ylim()[::-1])\n",
    "pylab.tight_layout()\n",
    "pylab.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/stewart/.virtualenv/PySeis/lib/python2.7/site-packages/ipykernel/__main__.py:17: RuntimeWarning: invalid value encountered in greater\n",
      "/home/stewart/.virtualenv/PySeis/local/lib/python2.7/site-packages/matplotlib/pyplot.py:423: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_num_figures`).\n",
      "  max_open_warning, RuntimeWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "777.01802206\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from matplotlib import collections\n",
    "\n",
    "def func(panel): \n",
    "    \n",
    "    panel['trace'][:,-1] = np.nan\n",
    "    trace_centers = np.linspace(1,284, panel.size).reshape(-1,1)\n",
    "    scalar = 284/(panel.size*0.5)\n",
    "    y = np.meshgrid(np.arange(1501), np.arange(284))[0].ravel() \n",
    "    offsets = (np.meshgrid(np.arange(1501), np.arange(284))[1]+1).ravel()\n",
    "    x = ((panel['trace']*scalar)+trace_centers).ravel()\n",
    "    \n",
    "    fig,ax = pylab.subplots()\n",
    "    #or i in range(284):\n",
    "          \n",
    "        #ax.plot(x[i],y[i],'k-')\n",
    "    ax.fill_betweenx(y,offsets,x,where=(x>offsets),color='k')\n",
    "\n",
    "    pylab.xlim([0,284])\n",
    "    pylab.ylim([0,1500])\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    pylab.tight_layout()\n",
    "    pylab.draw()\n",
    "    pylab.clf()\n",
    "\n",
    "dmap = np.memmap(file, dtype=toolbox.typeSU(1501), mode='r')\n",
    "eps = np.unique(dmap['ep'])\n",
    "for ep in eps[:1]:\n",
    "    panel = dmap[dmap['ep'] == ep].copy()\n",
    "    panel = toolbox.agc(panel, None, **params)\n",
    "    t = Timer(\"\"\"func(panel)\"\"\", setup=\"from __main__ import func; from __main__ import panel\")\n",
    "    print t.timeit(100)\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
